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Abstract 

The effective force between two parallel DNA molecules is calculated as a 
function of their mutual separation for different valencies of counter- and salt 
ions and different salt concentrations. Computer simulations of the primitive 
model are used and the shape of the DNA molecules is accurately modelled 
using different geometrical shapes. We find that multivalent ions induce a sig- 
nificant attraction between the DNA molecules whose strength can be tuned 
by the averaged valency of the ions. The physical origin of the attraction is 
traced back either to electrostatics or to entropic contributions. For multi- 
valent counter- and monovalent salt ions, we find a salt-induced stabilization 
effect: the force is first attractive but gets repulsive for increasing salt con- 
centration. Furthermore, we show that the multivalent-ion-induced attraction 

does not necessarily correlate with DNA overcharging. 
PACS: 87.15.Kg, 61.20.Ja, 82.70.Dd, 87.10.+e 
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I. INTRODUCTION 



During the last decade the question that concerns the possible existence of long-ranged 
attractive interactions between similarly charged objects in electrolyte solutions has been 
intensely debated. Experimental evidence of such an attraction is seen for deoxyribose 
nucleic acid (DNA) molecules [1-6], colloidal rods [7], charged clay particles [8], charged 
microspheres [9-11] and charged plates [12,13]. In particular, DNA molecules in solution 
are a paradigm for negatively charged polyelectrolytes due to ionization of its acidic phos- 
phate groups [14]. The DNA conformations display a considerable sensitivity to the ionic 
surrounding. The mutual repulsion of DNA polyions has to be overcome to form compact or 
condensed DNA bundles. Experiments show that DNA condensation occurs when about 90 
percent of its charge is neutralized by condensed counterions [1,3,15,16]. Such a strong neu- 
tralization of the DNA charge could be achieved by divalent and higher-valent counterions 
[17,18]. Besides of the phosphate neutralization, the multivalent ions induce an additional 
attraction between the DNA macroions mediated by strong correlation effects [16,19-26]. 
Thus the small ions play a complex role in DNA-DNA interactions and are not simply agents 
to screen the long-range electrostatic interaction. For example, they adsorb onto the DNA 
surface and can create bridges between the DNA molecules at small DNA-DNA separations, 
resulting in an ion cross-link attraction [24,27]. 

The electrostatic interaction between highly charged polyelectrolytes is usually treated 
within the framework of classical double- layer theory [28]. This theory is based on the 
mean-field Poisson-Boltzmann equation [29-36] and predicts a repulsion between similarly 
charged macromolecules. Though different modifications of Poisson-Boltzmann theory have 
been developed to account for ion-ion hard core correlations [37-40], an attractive contribu- 
tion in the double-layer theory is usually introduced via the van der Waals interaction forces 
[9,41-43]. However, the van der Waals forces alone cannot explain the experimentally ob- 
served attraction, since the Hamaker constant extracted from the experiments is artificially 
high [9,44-46]. 
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Theoretical investigations and numerical simulations indicate that an attraction between 
similarly charged objects emerges beyond the mean-field approaches. It is now a well estab- 
lished fact that the charge correlations and fluctuations in highly charged electrolytes can 
induce an attraction between the macroions [23,47-69]. Due to the very resemblance of the 
short DNA fragments to charged rods, the latter is a widely used toy model for the DNA 
molecule in theoretical treatments and computer simulations [24,30,63,70-75]. However, the 
details of DNA, such as the discreteness and helical structure of the DNA phosphate charges 
and the grooved shape of the DNA molecule, become essential as one approaches its sur- 
face. In this case, strictly speaking, all atom DNA simulations in molecular water would be 
a proper choice [76]. Unfortunately such sophisticated simulations can only be applied to 
small systems and small salt concentrations [77]. Thus, first, a devision of a "sophisticated" 
DNA model, which goes beyond the simple homogeneously charged cylinder model, and 
second, an investigation of the interaction forces between such DNA molecules, remains a 
challenging task. 

This paper is an extension of our previous works on DNA electrostatics [78,79]. In 
Ref . [78] , simulation results for the DNA-DNA interaction were compared with the predic- 
tions of different linear theories for the case of monovalent ions. It was shown that the 
DNA-DNA interaction, at separation distances smaller than the Debye screening length, 
differs from the predictions of mean-field theories. This provides evidence that the inter- 
molecular interaction depends not only on how many ions are in the DNA proximity (which 
is exactly what the ordinary linear theories rely on), but also on where those ions are lo- 
cated relative to the DNA structure, i.e., whether they penetrate into the grooves or not. 
In Ref. [79], on the other hand, a detailed distribution of ions of different valencies and mo- 
larities near the DNA surface was explored for a more realistic, grooved shape of the DNA 
molecule. The results obtained indicate that the paths of counterion and colon condensa- 
tions strongly depend on the DNA surface geometry. Taking this into account we expect 
that the implemented DNA models with different geometries will also affect significantly 
the effective DNA-DNA interaction. Thus, in this paper, we will focus on the mechanism 
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of attraction between two DNA molecules shaped similar to models introduced in Ref . [79] . 
Our goal is to see the effects, which increasing detail of various DNA models have on the 
DNA-DNA interaction. We show that the DNA shape is an essential contributor to the in- 
teraction force for multivalent counterions, whereas it has a minor effect on the interaction 
force for added multivalent salt. The origin of the attraction in the simple and sophisticated 
DNA models is different. For instance, a Coulomb depletion-like attraction [80] for the 
salt-free case depends on the implemented DNA model. It has been revealed that there is a 
non-monotonic force-salt dependence at a fixed DNA-DNA separation for added monovalent 
salt and divalent counterions. This is exemplified by the variation of the interaction force 
from a strong attraction towards a strong repulsion and a following decrease in magnitude. 
Detailed investigations connect this "salt-induced stabihzation" to the entropic part of the 
total interaction force. We also address the competition between the multivalent counterion 
and the multivalent salt-induced attractions. It is shown that the increase of the divalent 
salt concentration at a fixed monovalent ion number drives the DNA-DNA interaction force 
into an attraction through the overcharging of DNA molecules. However, the DNA-DNA 
attraction induced by trivalent counterions decreases, while the DNA molecule gradually 
gets overcharged due to added divalent salt. 

The reminder of paper is organized as follows. We give a short general overview of ion 
binding and DNA condensation in Section II. The system parameters and quantities studied 
in the present work are discussed in Section III. Sections IV and V contain simulation 
details and the implemented simulation techniques. The specific DNA configurations at 
short DNA-DNA separations are discussed in Section VI. Sections VII and VIII are devoted 
to simulation results for monovalent and divalent salt ions respectively. We conclude in 
Section IX. 
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II. ION BINDING AND DNA CONDENSATION 



There are essentially two contenders for the dominant attractive force in the DNA con- 
densation: hydration forces [81,82] and correlated counterion charge fluctuations. Through- 
out this paper we neglect the granular nature of water and the solvent- induced forces [83] , 
and concentrate only on the electrostatics of the DNA condensation. The water dielectric 
effects and hydration forces will be briefly (and quahtatively) discussed in Section IX. 

Under physiological conditions, the DNA molecule is surrounded by an ionic atmosphere 
with a Debye screening length in the range of SA-IOA. Within the distances r < Xd 
above the DNA surface, a nonhnear screening of the DNA phosphate charges takes place. 
Hence, if the surface-to-surface separation between two DNA molecules is less than Xd, a 
nonlinear theory [84-87] has to be apphed. At surface-to-surface separation distances on 
the order of or beyond A^, Debye-Hiickel theory (based on a linearized Poisson-Boltzmann 
treatment) is a reasonable approximation to describe the ionic atmosphere around the DNA 
molecule. In Ref. [78] we have examined several mean-field theories for their ability to match 
the numerically calculated DNA-DNA interaction forces: the homogeneously charged rod 
model, the Yukawa segment (YS) model, and the Kornyshev-Leikin (KL) theory [88]. For 
the case of an overall monovalency of counterions and salt ions, both the simulations [78] 
and the above mentioned theories reveal repulsive forces between the DNA molecules for 
all mutual orientations and separation distances. We have shown that, except for short 
separation distances, there is a qualitative agreement between the theoretical and numerical 
results if a proper charge and size renormalization in the former is performed. 

For multivalent counterions and added salt ions, there is experimental evidence that the 
DNA molecules attract each other. Such an attraction is completely missed in the hnear 
theories such as the homogeneously charged cylinder and YS model. In these theories all the 
nonlinear salt effects are again accounted for through the phosphate charge and screening 
length renormalization procedure. Only the mean-field KL theory predicts a DNA-DNA 
attraction for some DNA-DNA separations and azimuthal molecular orientations. In detail. 
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the KL theory distinguishes between strongly condensed (also called as bound or adsorbed) 
and a cloud of diffusive (non-bonded) counterions. A tight adsorption is assumed to take 
place within the Stern layer of thickness ^ = 74/47rAs=2A, where = e^/{ekBT) is the 
Bjerrum length, A is an average area per elementary charge on the DNA surface, e is a 
dielectric constant of solution and ksT is the thermal energy. The KL theory [88] predicts 
an attractive force between the two DNA molecules if the following conditions are fulfilled: 
i) more in-groove than on-strand condensation, ii) the right complementary alignment of 
the positively charged grooves on one helix facing the negatively charged strand on the 
other helix. In other words, the KL theory assumes that it is the DNA charge helicity 
that entails an intermolecular attraction for surface-to-surface distances in the range of 
8 — ISA. Theoretical results and computer simulations [35,78,89-94], however, indicate that 
no charge-helicity effects extend further than few A from the DNA surface. There is also 
experimental evidence [95] that at surface-to-surface separation distances comparable with 
the Debye screening length, the DNA-DNA separation does not affect the DNA orientation. 

The discreteness of the DNA phosphates, explicitly taken into account by our DNA 
models, enhances the counterion concentration [96] and the surface adsorption of ions [97] 
through the increased Coulomb coupling between the phosphates and the counterions. This 
boosts the counterion correlations near the DNA surface [98]. Experiments indicate that 
the divalent counterions, depending on their in-groove or on-strand localization [99], have 
different impact on the DNA systems. Thus, the transition metals with higher affinity to the 
DNA bases [100,101] condense on DNA [102,103], while alkah metals do not [4,103]. On the 
other hand, the chemical identity of the cation is a factor of minor importance compared with 
the magnitude of their charge when gc > 2 [3,4,16,104-107]. Thus the spermidine Spd3+ and 
spermine Spm4_|_ ions, abundant in living cells [108-111], neutralize the negatively charged 
DNA backbone predominantly via the non-specific (Coulomb) interaction [112-114]. This is 
supported by new experiments [21,115], polyelectrolyte and counterion condensation theories 
[116-119], and computer simulations [89,120]. 
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III. SYSTEM PARAMETERS 
A. DNA Models. 

The B form of DNA has an inner core of radius 9A formed by nucleotide pairs, and two 
sugar-phosphate strands spiralUng around it. The latter form the well-known double helix 
with a pitch length P about 34A [14]. There are two phosphate groups per base pair, and 
10 base pairs per pitch length, or hehcal turn. The axial rise per base pair in the DNA long 
axis is 3.4A, thus there is one elementary charge per each 1.7A [121]. The average value 
of the angle between the adjacent base pairs is 36° and the average distance between the 
neighboring charges on the DNA surface is about 7A. This distance is much smaller than the 
helical pitch and of the order of Debye screening length under the physiological conditions. 
There is a small shift in the z coordinate of two opposing phosphates belonging to different 
hehces of DNA, 5z = 0.34A. 

Three DNA models, a cylinder model (CM), an extended cylinder model (ECM) and 
the Montoro-Abascal model (MAM), are considered. Our aim is to obtain a detailed under- 
standing of the physical mechanism of ion-mediated DNA interactions, in particular how the 
geometry of different DNA models gives rise to new effects. The CM has a hard cylindrical 
core of diameter D — 20A and two strings of monovalent phosphates of size dp = 0.4A. 
The KL theory, and almost all the Poisson-Boltzmann like theories and most of primitive 
model (PM) computer simulations, have utilized the CM as a simple DNA model. In the 
ECM, first designed by Lyubartsev et. al. [89], the helical grooves of DNA are incorporated 
through the shrinking of the DNA core to the size D — 17.8A and sweUing the phosphate 
spheres to the size dp = 4.2 A. A grooved structure, which resembles the real DNA appear- 
ance, is achieved in the MAM [90] through the adding another neutral sphere between the 
cyhndrical core and the charged phosphate sphere. The cyhndrical core in the MAM has 
a diameter D = 7.8A, the inner string of neutral spheres is centered at a radial distance 
r = 5.9A, and the outer string of phosphates is centered at a radial distance r = 8.9A. 
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Both spheres have the same and z coordinates and diameter dp = 4.2A. A full description 
of these models is given in Refs. [78,79,90]. 

In addition to the two DNA molecules the system contains counterions of charge 
symmetric salt ions of concentration Cg and charges and All the small ions are 
modelled as a hard spheres of a diameter dc for counterions, c?+ and d- for the salt ions. 
The whole system is held at room temperature T = 298i^. The primitive model simulations 
with no explicit water deal only with a passive (non-specific) binding and completely neglect 
the specific binding of counterions to the DNA grooves. In this case the ion binding sites 
are determined by the steric and Coulombic interactions [79,122]. 

The interactions between the mobile ions and the phosphate charges are described within 
the framework of primitive model as a combination of the excluded volume and Coulomb 
interactions reduced by the inverse of the dielectric constant e of the solvent. The corre- 
sponding pair interaction potential between the different charged hard spheres is 

foo for r < (di + dj)/2 

ViAR) = \ , . (1) 

\mi_ for R> {d, + d,)/2 
where R is an interparticle separation distance, i and j are indices denoting the different 
particles species. Possible values for i and j are c (for counterions), -|-, — (for positively and 
negatively charged salt ions) , p (for phosphate groups) and n (for the neutral spheres in the 
MAM with gji=0). In addition, there is an interaction potential between the DNA hard 
cylinder and the free ions i — c, -|-, — . This potential has a simple excluded volume form 
such that the free ions cannot penetrate into the cylinder. 

B. Simulated Quantities. 

Our basic simulated quantity is the effective force [78,80] between the DNA molecules 

F = Fi + F2 + F3. (2) 

Here Fi is the direct Coulomb force acting onto all the phosphate charges belonging to one 
helical turn of one DNA molecule as exerted from the phosphate groups of the other DNA, 
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^i = -E V,-. E V^{\r^,-rZ\)\. (3) 

k \ n=l;n^k ) 

The sum Y!^. only runs over phosphates of one hehcal turn of the DNA molecule. 

— * 

The second term F2 corresponds to the Coulomb interactions between the phosphate 
charges and the mobile salt ions. This term describes the screening of the DNA charge, 

F2 = -E ((EE ^r^M\ - D) ) ■ (4) 

k \ j=c,+,— ;=i / 

The third term F3 arises from the entropic contribution of small ions due to their excluded 
volume interaction with the DNA molecular surface Si. Its value for one helical turn is 



Fs^-ksTf df I E 



pA^ , (5) 



where / is a surface normal vector pointing outwards the DNA cylindrical core. This term 
becomes increasingly important as the Coulomb coupling parameter Fpc is elevated for the 
multivalent counterions [78,80,123,124], 

= I f I , , , ■ (6) 

Qc dp + dc 

The parameter Tp^ determines the importance of thermal fluctuations. When Tp^ > 1, the 
Coulomb interaction energy between the DNA and the surrounding salt ions dominates over 
the thermal fluctuations in system. 



IV. COMPUTER SIMULATIONS 

We consider two parallel DNA molecules, separated by a distance R along the xy diagonal 
of the cubic simulation box of size L and volume V — L"^. The size of the simulation box 
L = IO2A corresponds to the three full turns of B-DNA [90]. The box also consists Nc 
counterions and — — Ng salt ions of both signs. The counterion concentration 
is fixed by the charge of DNA molecules in the simulation box due to the constraint of 
global charge neutrality. A typical snapshot of the simulation cell is illustrated in Fig. 1. 
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Periodic boundary conditions in all three directions are applied to reduce the confined volume 
effects in electrolytes. The DNA replicas in z direction produce an infinitely long DNA 
molecule which avoids the end effects encountered in other molecular simulations of short 
DNA segments [125-127]. The phosphate spheres are monovalent, qp — — e, where e is 
an elementary charge. The total number of simulated salt ions is varied from to 2000 
depending on the salt concentration in the bulk. Hereafter we will refer to each simulation 
by its nominal salt concentration Cg defined as a ratio between the total number of ions Ng 
in the cell without the DNA molecules and the system volume V ^ Cg = Ng/V. The dielectric 
permeability e is considered to be a constant everywhere in system, which avoids the need of 
electrostatic images [128,129]. The long range interactions between the two charged species 
and their replicas in the neighboring cells are handled via the Lekner method [130] and its 
modification for the particular cases, when pair of interacting charges are sitting exactly on 
one of the coordinate axies (this case was considered in Appendix of Ref. [78]). 

We have performed extensive molecular dynamics (MD) simulations for a range of dif- 
ferent microion valencies. The simulated states are given in Table 1. The ion diameter was 
chosen to be dc=3A. This parameter defines the closest approach of the ion to the DNA sur- 
face and has a strong impact on the polyion electrostatics [20,131,132]. A test simulation for 
an increased ion diameter dc=hK, which mimics the ion hydration in solvent [20,79,93,132], 
shows no qualitative changes of the reported results. We want to emphasize again that 
specific ion effects, as exemplified by the Hofmeister effect [133,134], are not accounted for 
in our model. 

During the simulations, we calculate the interaction forces between the two DNA 
molecules for different separation distances R. Due to the strong screening of the DNA 
phosphates, the actual salt concentration in the bulk of the simulation box C^(i?), measured 
far away from the polyelectrolytes, is R dependent and is always smaller than the nominal 
salt concentration Cg, C'g{R) < Cg. Thus, an implementation of the conventional MD pro- 
cedure with a fixed ion number Ng will yield to the interaction forces which correspond to 
bulk densities C'g{R) for each intermolecular distance R. This problem can be avoided by 
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considering a solution with a constant a chemical potential // via the Grand Canonical (GC) 
simulation method [67]. The GC simulation is a natural choice to mimic the experimental 
situation where the actual salt concentration of the ordered DNA phase is not known a 
priori. Instead, it is given by the thermodynamic condition that the chemical potential 
in the DNA solution has to be the same as in the bulk electrolyte phase with which it is in 
equilibrium and whose concentration Cg is experimentally known. Thus the number of ions 
in the simulation cell is automatically adjusted to the specified value of chemical potential 
/X, which, in turn, is linked to the concentration of ions in the bulk phase Cg [135-138]. In 
the present paper a combination of different grand canonical molecular dynamics (GCMD) 
schemes is used which is optimally suited for our task. 

V. GRAND CANONICAL MOLECULAR DYNAMICS 

In addition to the usual propagation of the particles, the conventional GC simulation 
technique [137,138] consists of the creation of particle at a random position in the simulation 
box or destruction of a randomly chosen particle. Each of these moves is associated with 
a probability of acceptance, which is determined by the ratio of two Boltzmann factors. 
In application to electrolytes a modified GC method was devised in Refs. [71,120,139-142], 
where the insertion or removal of a pair of ions of the same valency and opposite charge 
is done simultaneously to keep the system electroneutral. Unfortunately, these moves have 
relatively low acceptance rates for a dense and multivalent salt solutions [71], making the 
simulations inefficient. Special methods, similar to the cavity-biased method [141] and the 
gradual particle insertion method [143] developed for uncharged systems, have to be applied 
to overcome these obstacles in electrolytes. A biased insertion/destruction procedure, apt 
for an application to low temperature ionic fluids, is reported in Ref. [144]. 

Another challenge in the GC simulations is the apparent incompatibility of the deter- 
ministic and stochastic approaches. The dynamical information will be adversely affected 
when particles suddenly appear and disappear. This effect becomes even more pronounced 
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for a non-homogeneous system [55,56], like a DNA immersed in solution, where an artificial 
and unreahstic ion flow toward the DNA surface appears. This will further destabilize the 
system equilibrium. To minimize this inconsistency of the system dynamics, a method of 
local potential control (LPC), first introduced in Ref. [145], can be adopted. Within the 
LPC method, the creation and destruction of particles is restricted to a control volume. 
The other possibility is a procedure developed by Attard in Ref. [146] , where the GCMD is 
performed with a fixed number of particles by coupling the variations in the system size to 
the instantaneous chemical potential determined by the virtual test particle method. This 
method also cures the low acceptance rates of particle insertions and deletions for dense 
systems. In the present simulations we take advantage of both the above mentioned meth- 
ods [145,146]. In detail, we first determine the specified nominal chemical potential /j, of 
the bulk electrolyte in the absence of DNA molecules via a modified Widom method with 
multiparticle insertion [147,148]. Then we match the actual chemical potential to the 
nominal chemical potential n using a GCMD simulation similar to the method invoked in 
Ref. [145] and locate the control volume near the cell boundary. At each time step an equal 
but arbitrary number of creation/deletion attempts are made in the control volume. Af- 
ter a successful creation, a velocity is drawn from the Maxwell-Boltzmann distribution at 
a temperature T and assigned to the new particle. In the first stage the particle number 
Ng in the simulation box increases monotonically from its initial value Ng given in Table I. 
Then N'^ approaches its saturated value and starts to fluctuate around it. This is followed 
by the fluctuation of instantaneous chemical potential /x' around the /x. At this stage we flx 
the particle number and allow the system size to fluctuate according to procedure given by 
Attard in Ref. [146]. The fluctuations along the x and y directions {z direction is strictly 
bound to the DNA length) never exceed a few percents of the box size L. Our test simu- 
lations with and without the Attard method [146] show the equivalence of the algorithms, 
with the former being much faster. 
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VI. MUTUAL DNA CONFIGURATIONS 



We calculate the total interaction force F{R), Eq.(2), and its components F2{R) and 
F^lR), compare Eqs. (4) and (5), respectively, for a given nominal salt concentrations 
Cs- The direct phosphate- phosphate interaction Fi{R), Eq.(3), does not depend on the 
salt density and its assessment is straightforward. It should be mentioned that in addition 
to the separation distance R, there are three angular variables which define the mutual 
configuration of two DNA molecules. These variables, the azimuthal angles 4>s, 4>o and 
are shown in Fig. 2. The angle 0^ defines the widths of the DNA grooves [79] in the xy 
plane, it is 134° for the CM and ECM, and 154° for the MAM. The parameter 0o is the 
angle between the phosphate charge and the DNA-DNA separation vector R = Ri — R2 
and characterizes the discrete location of the phosphate charges along the strands. All the 
results for the DNA-DNA interaction are periodic in (f)Q with a periodicity of 36°. The angle 
4> describes the rotation of the second DNA cylinder around its long axis with regard to 
the first DNA cylinder. There are five particular DNA-DNA configurations which make 
a strong contribution to the interaction force at short separation distances [78]. Three of 
these configurations correspond to the case when the phosphate charges of neighboring DNA 
molecules are "touching", see Fig. 3. The other two particular configurations correspond 
to the so called "DNA zipper" situation, when the strands of one DNA stand against the 
grooves of the neighboring DNA [149]. This happens when = ±37r/5 regardless the value 
of 00- Our previous [78] and present simulations prove that the interaction force does depend 
on the mutual DNA configurations at a short separation distances R < 25A, or when the 
surface-to-surface distance between the two DNA molecules is less than SA. This follows 
mainly from ion bridging between the two neighboring phosphates via a positive salt ion in 
configurations pictured in Fig. 3, or from sharing an adsorbed salt ion in one DNA groove 
and on the other DNA strand in the DNA zipper configurations. A short range attractive 
interaction between the charged rods arising from interlocking counterions, also sometimes 
called counterion cross-hnks, between the rods has been investigated in Ref. [24,150,151]. In 
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fact, for such small separations, comparable with the solvent particle size, discrete solvent 
effects will show up in vitro [82,97]. On the other hand, the multivalent ions increase 
the hydrodynamic radius of the DNA molecule, which in turn makes it unlikely for two 
neighboring DNA molecules to come closer than the contact shell-to-shell distance [132]. 
Arguments against the existence of cross-links for multivalent ions are given in Ref. [1] 
in order to explain the fluidity of the condensed phase of DNA system. X-ray scattering 
experiments in DNA aggregates a show that the surface-to-surface spacing between the DNA 
helices is only about one or two water molecule diameters [150]. Thus, numerical results 
for small DNA separations, accessible in simulations but subject to a complicate statistical 
averaging procedure, bear no physical meaning to match the experimental results. For larger 
separation distance, R > 25A we find no detectable dependence of the interaction forces on 
the azimuthal angles 4>q and 4>. This is in accordance with the early reports [35,89-94] that 
the hehcity and discreteness effects of the DNA charges are generally small and dwindles a 
few angstroms away from the DNA surface. In all figures hereafter we show the interaction 
forces starting from the distances R = 24A. The interaction forces are scaled per DNA 
pitch, i.e. per 10 DNA base pairs. A positive sign of the forces denotes a repulsion, while a 
negative sign denotes an attraction. The cases of monovalent and multivalent salt ions are 
considered separately. 



VII. SIMULATION RESULTS FOR MONOVALENT SALT. 

A. Monovalent Counterions. 

The calculated DNA-DNA interaction forces for monovalent salt and counterions are 
depicted in Fig. 4. All three DNA models exhibit a repulsion between the DNA molecules 
for all calculated separation distances and salt densities shown in Fig. 4 [2,78]. The repulsion 
in the CM is roughly twice as strong as in the MAM. This is a result of grooved nature of 
the MAM where the vast majority of adsorbed counterions sits in the minor groove [79]. 



14 



The salt dependence of the force at a fixed, small separation appears to be non-monotonic. 
This is clearly visible only for the ECM and MAM but hardly detectable for the CM. In 
detail, if the salt molarity is increased from Cs=0 mol/1 to Cs= 0.024 mol/1, the repulsion at 
short distances becomes stronger, see the inset of Fig. 4c. Though this trend is at odds with 
the classical screening theories, a similar effect has already been reported in Ref. [78,124]. 
A detailed consideration of the interaction force components in Fig. 5 reveals that the 
non-monotonicity of the interaction force F{R) has a purely electrostatic origin. Both the 
electrostatic force F2{R) and the entropic force Fs{R) are repulsive for all indicated salt 
densities; the non-monotonicity is only contained in F2. The non-monotonicity is a tiny 
effect for monovalent ions, but for the divalent counterions and monovalent salt, this non- 
monotonicity aggravates and induces a switch from attraction to repulsion, see the next 
subsection. 



B. Divalent Counterions. 

The interaction forces for divalent counterions and monovalent salt, given in Fig. 6, reveal 
a DNA-DNA attraction for small added salt concentrations. First we analyze the salt free 
case, when the attraction in the CM is nearly three times stronger than in the MAM. The 
origin of these attractions in the CM and MAM is completely different. In the cylindrical 
model the attraction is totally associated with the "Coulomb depletion" effect [80,152]. Such 
an ion depletion is related to the formation of strongly correlated counterion liquid on the 
DNA surface. According to the results of Ref. [79] , the dense spots of counterion hquid are 
mostly associated with the DNA phosphates. For short DNA-DNA separations, when the 
mean separation distance between counterions on the DNA surface exceeds the separation 
distance R, the two strongly correlated counterion clouds on the different DNA rods start 
to repel each other. As a result, the local density of ions in the inter DNA area becomes less 
than its value on the outside area. This in turn leads to a unbalanced pressure from the outer 
counterions [153], which implies an attraction, as proven by the solid lines in Fig. 7. It is 
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worth to mention that in Ref . [80] such a correlation mechanism was reported for a spherical 
colloid with a central charge in a low dielectric medium. Surprisingly, for the cylindrical 
macroions with a discrete surface charge considered here, we recapture a similar effect. Note 
that our results do not support the implications of Ref. [118], where an attractive force 
between the two rodlike polyions is assumed to be mediated by the sharing of condensed 
counterions. 

Contrary to the CM, the attraction in the MAM and ECM has a purely electrostatic 
origin, as proven by Fig. 8. There are more ions in the inter DNA-DNA area compared to the 
outer DNA-DNA area. The range and strength of this attraction is higher for the ECM than 
for the MAM. To our believe, this is due to the different counterion condensation patterns 
on the DNA surface in the MAM and ECM. In the latter model the ions predominantly 
adsorb to the DNA strands and in the minor groove. Therefore they occupy more DNA 
surface area compared to the MAM, where the main destination of the ion adsorption is the 
minor groove of DNA [79]. According to a simple intuitive picture, the adsorbed divalent 
counterions form a strongly correlated fluid on each DNA surface. A mutual arrangement 
of these two neighboring shells with a minimal potential energy gives rise to the attractive 
electrostatic force between the DNA molecules. 

Upon an addition of monovalent salt, the trend in the interaction force depends sensi- 
tively on the DNA shape which is modelled differently in the CM, ECM and MAM. In the 
CM the DNA-DNA attraction persists to high salt concentrations. Only at short separations 
and dense salt the interaction has a repulsive branch. However, there is a counterintuitive 
behavior of the force-salt dependence in the ECM and MAM: a small increase in the salt 
concentration Cs suppresses the attractive interaction force F. As Cg increases further, the 
DNA-DNA interaction force becomes strongly repulsive over a broad range of separation 
distances. At even higher C^, the interaction force F is completely screened out and de- 
scends toward zero in accordance with the classical double layer theories. We denote this 
trend salt-induced stabilization. It is in complete contrast to salt-induced dcstabilization 
or salt-induced coagulation which is typical for charged colloids [154-156]. The force-salt 
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dependence at a fixed separation distance is shown in the inset of the Fig. 6c. To explain the 
mechanism of salt-induced stabilization, we separately plot the interaction force components 
in Fig. 7 and Fig. 8 for the CM and MAM, respectively. As the salt density is increased, 
the entropic force in the CM first changes its sign from an attraction to a repulsion, then 
falls back to zero from above. The electrostatic force in the CM behaves in a similar manner 
but with opposite sign. It first goes down from the positive to the negative values, then 
approaches zero from below. Hence, at higher salt densities the attraction, seen in Fig. 7, is 
electrostatically driven. 

Contrary to the CM, the non-monotonicity of the interaction force in the ECM and MAM 
has an entropic origin, see Fig. 6 and Fig. 8. For a fixed separation distance it! it first goes up 
and then drops back to zero as Cg is increased. Similar to the monovalent counterion case, 
the entropic force is repulsive over the full range of separation distances. The electrostatic 
force, which was totally repulsive for the monovalent counterions, now is attractive for dilute 
salt densities. Thus, it is the electrostatics that makes the total interaction attractive for 
the ECM and MAM at small salt densities. As Cg is increased, the monovalent ions tend to 
replace the adsorbed divalent ions on the DNA surface in accordance with the two- variable 
theory of Manning [157]. This replacement is energetically favorable, since the divalent ions 
gain more polarization energy in the bulk electrolyte. This results in the loss of the attractive 
electrostatic force and the weakening of the entropic force. At the final stage, when all the 
divalent sites on the DNA surface are occupied by monovalent ions, the strongly correlated 
fluid structure is destroyed and the entropic and electrostatic forces drop to zero. 

An other interesting observation is the occurrence of zero force in the CM in Fig. 6 at 
small separation distances. As the salt concentration is increased, the distance where the 
interaction force vanishes, shifts towards larger values. The physical meaning of this effects 
is not clear to us. Probably it is an artifact of the simple DNA model with no grooves. 

We note that the salt densities invoked in the current simulations are below the con- 
centrations where a salt depiction [155] appears near the DNA surface. Thus the observed 
attractive DNA-DNA force is not induced by the salt depletion effect, which was claimed in 
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Ref. [155] to be one of the main contributors to the protein-protein attraction. Furthermore 
the attraction observed here is unrelated to metastable ionization [98]. 

C. Trivalent Counterions. 

Overwhelmingly attractive DNA-DNA interaction forces for trivalent counterions and 
monovalent salt ions are plotted in Fig. 9. This attraction is stronger than for the divalent 
counterions, in accordance with the results of Ref. [63]. Similar to the divalent counterion 
case, the DNA-DNA attraction in the CM is much stronger than the attraction in the MAM 
for a given salt density. Evidently this is related to the different counterion condensation 
patterns for different DNA models [79]. There is no salt-induced stabilization for the ECM 
and MAM for the salt densities indicated in Fig. 9. Test runs with higher salt concentrations 
(not shown here) reveal that the salt-induced stabihzation in the MAM does appear at 
Cg — 1.2 mol/1 when the repulsive interaction force starts to descend towards zero. However, 
up to these high salt concentrations, the attractive force in the ECM steadily approaches 
zero from below. The same trend holds in the CM as well. 

The entropic force is non-monotonic against added salt both in the ECM and MAM, 
see Fig. lib. The reason is the following. The added salt enhances the release of condensed 
trivalent counterions from the DNA surface to the bulk [62]. However, the strong electro- 
static correlations in the inner area between the DNA molecules hinder this release. As a 
result, the trivalent counterion release will be asymmetric over the surface of DNA molecule. 
This will lead to an excess entropic force that pushes the DNA molecules away from each 
other. The osmotic pressure of added salt balances and at some salt density overcomes this 
entropic force, eventually driving the entropic force to zero. This entropic non-monotonicity 
does not survive in the total interaction force shown in Fig. 9. This is because of the strong 
electrostatic attraction between the DNA molecules. Whereas the entropic repulsion is of 
the same order for divalent and trivalent counterions, the electrostatic attraction for trivalent 
counterions is roughly twice as strong as for divalent counterions. 
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The attraction in the CM for small added salt densities Cg in Fig. 10 has again an entropic 
origin [80]. For higher salt densities, however, the attraction in the total interaction force 
at the short separation distances is electrostatically driven. 

VIII. SIMULATION RESULTS FOR MULTIVALENT SALT 
A. Monovalent Counterions. 

The results obtained in previous sections indicate that the multivalent counterions gen- 
erate strong correlations inside the system and induce an electrostatic attraction between 
the DNA rods in the monovalent salt system. An intriguing question is how this DNA-DNA 
attraction relates to DNA overcharging. Going back to the single DNA case, considered in 
Ref. [79], we observed no overcharging for multivalent counterions and monovalent salt ions. 
Thus we conclude that the overcharging effect is not a necessary condition for a DNA-DNA 
attraction to take place. In other words, the electrostatic ion correlations, which are not 
strong enough to induce the macroion overcharging, are able to induce an attractive inter- 
molecular force. On the other hand, in Ref. [79] we have seen a DNA overcharging when 
multivalent salt ions were pumped into the DNA suspension. When a DNA overcharging 
occurs, the ions form a sequence of radial alternating charged layers around the DNA sur- 
face. The width of these layers is comparable with the Debye screening length Ad. The 
question we want to address here is to what extent the existence of such layers affects the 
DNA-DNA interaction. 

The total interaction forces for the divalent salt and monovalent counterions are depicted 
in Fig. 12. For all three DNA models, the DNA-DNA interaction is repulsive at small 
salt densities Cg and transforms to an attraction for a sufficiently high C^. The DNA- 
DNA repulsion at lower salt is composed from the repulsive F2 and F3 components of the 
interaction force F. In a similar way the attraction at a dense salt, which is strong for 
the CM and weaker for MAM, arises from both attractive electrostatic and entropic forces. 
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For like-charged colloid particles at similar parameters for the small ions, an attraction was 
reported only for the electrostatic component of the interaction force [158]. A multivalent 
salt-induced precipitation of polyelectrolyte solution is also addressed in Ref. [159]. Whereas 
the van der Waals or hydrophobic forces could be presumed as a source for this attraction 
[46,81,160], our simulations are able to catch this effect without resorting to these forces. 

The DNA-DNA attraction at the divalent salt density Cs=0.71mol/1 in Fig. 12 corre- 
sponds to the overcharging of a single DNA molecule, see Ref. [79] . Thus, an overcharging 
and entailed charged layers near the DNA surface correlate with the DNA-DNA attraction 
in the divalent salt, at least when monovalent counterions are involved. This conclusion con- 
tradicts the results of Refs. [161] where the onset of the DNA overcharging was considered 
to entail a repulsion between the DNA molecules, and therefore, a reentrance of the DNA 
condensation. From our point of view, the discrepancy between our result and the result of 
Ref. [161] is due to the different definition of overcharging. Whereas we count all charges 
near the DNA surface, only big multivalent ions were counted in Ref. [161]. 

B. Trivalent Counterions. 

Simulation results for divalent salt and trivalent counterions are illustrated in Fig. 13 for 
the MAM. Now the DNA-DNA interaction force and both of its components F2 and F3 (not 
shown here) are strongly attractive for all the calculated salt densities. The CM and ECM 
models exhibit similar trends. We note that at small salt densities, where no overcharging 
was found for a single DNA molecule [79] , the obtained attractive force relates to the strong 
charge correlations in system. Thus the claim that the overcharging effect is a sole factor 
that induces the DNA attraction, is not correct. Broadly speaking, there is a competition 
between the correlations that induce an attraction between the DNA molecules (multivalent 
counterion induced correlations) and the overcharging effect, which induce a DNA-DNA 
attraction (multivalent salt-induced correlations) . To understand the physics of this compe- 
tition we have analyzed the tendencies of these two correlation effects against the increase 
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of the salt density. Fig. 13 shows that the interaction forces decay monotonically with in- 
creasing salt concentration C^. This trend is in contrast to the results for a divalent salt and 
monovalent counterions in Fig. 12, where more salt- induced more attraction. Thus, the main 
contribution to the DNA-DNA attraction comes from the strong correlations between the 
two strongly correlated layers of trivalent counterions on the DNA surfaces. Pumping more 
divalent salt into the system destroys the two-dimensional crystal structure and correlations. 
But the interaction force remains attractive, mainly due to the additional overcharging of 
the DNA. As a result, the DNA-DNA attraction survives for a dense salt, in opposite to the 
case of trivalent counterions and monovalent salt shown in Fig. 9. 

IX. CONCLUSION 

We have studied the interaction forces between a pair of DNA molecules in an elec- 
trolyte that contains a mixture of monovalent and multivalent ions. Three models for the 
DNA shape, employed in our simulations, indicate the importance of the DNA geometry 
on the electrostatic and entropic forces in the DNA conformations. We show that the 
DNA-DNA attraction is related to the charge correlations in strongly charged systems. We 
distinguish between multivalent counterion and multivalent salt induced attractions. In gen- 
eral, the higher the mean valency of the microions in the solution, the stronger is the mutual 
attraction. Below we shortly summarize the main results of this manuscript. 

For the multivalent counterions, the DNA shape is an essential contributor to the inter- 
action forces. Thus: 

i) For no added salt the DNA-DNA attraction in the CM is related to the Coulomb depletion 
mechanism. This depletion effect results in an attractive entropic force. However such an 
attraction mechanism does not exist in the MAM. 

ii) For the nonzero added salt cases, an attraction in the CM is mainly due to a combination 
of electrostatic and entropic forces. However the attractive force in the MAM always has an 
electrostatic origin. The entropic force in the MAM is always repulsive. 
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iii) There is a non-monotonic force-salt dependence at a fixed separation distance. For diva- 
lent counterions, there is a change of the interaction force from the repulsion to an attraction, 
and then back to zero, which we call salt-induced stabilization. 

iv) An increase of the salt concentration suppresses the charge correlations and thus effec- 
tively screens the DNA-DNA attraction. 

For the multivalent added salt, the DNA shape has a minor effect on the interaction 
forces. DNA-DNA interaction forces are stronger for the CM than for the MAM. Further 
trends are: 

i) An increase of the divalent salt density at a fixed monovalent ion number drives the 
DNA-DNA interaction force towards attraction. Both DNA molecules are overcharged in 
the attractive force regime. 

ii) For trivalent counterions the addition of divalent salt decreases the DNA-DNA attrac- 
tion. The more the DNA becomes overcharged, the less is the attraction between the DNA 
molecules. 

iii) The correlation effects related to multivalent counterions have a greater infiuence on the 
DNA attraction than the correlation effects related to multivalent salt ions. In other words, 
an overcharging-induced attraction is weak compared to the counterion-induced attraction. 

We would like to make some comments about the range of the DNA-DNA attraction 
which directly influences the phase diagram of DNA solutions [162,163]. Compared to the 
Debye screening length, the attraction forces between the DNA molecules are long-ranged 
(they are beyond the screening length). In contrast, for colloids, usually the attraction is 
short ranged in comparison with the screening length [11]. Thus, the calculated attrac- 
tive DNA-DNA forces can lead to phase separation in DNA solutions but they are not the 
dominant driving forces of DNA condensation: The actual "arbiters" of phase instability 
in macroion solutions seem to be many-body interactions [164-166]. These forces separate 
the DNA assemble into DNA rich and DNA poor regions even for a purely repulsive pair- 
wise interaction [55,167,168]. Another intriguing effect is the appearance of a cholesteric 
phase of DNA condensates, which arises from DNA chirality and many-body interactions 
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[16,169-171]. 

Finally, let us comment on some limitations of our model. Our choice of a continuum 
dielectric model has the intention to separate the purely electrostatic effects from the effects 
of hydration and the molecular structure of the solvent. In many applications, including 
strong polyelectrolytes and high salt concentrations, continuum dielectric models (the prim- 
itive electrolyte model) have been successful (see Refs. [172,173]). However, strictly speaking, 
the continuum model is not justified at small ion-DNA and ion-ion separations where the 
molecular nature of the solvent is no longer negligible [83,174-176]. At these distances the 
effective (mean-force) potentials of ions have one to two oscillations around the potential of 
primitive model [91,97,175-178]. 

It is worth to mention the complete neglect of the specific binding (or chemisorbtion) 
of counterions to the DNA grooves in our simulations. Active binding [179] can be taken 
into account through the incorporation of a full water description [120,180-182] or via the 
implementation an additional sticky potential to certain ions in parts of the DNA areas. A 
cylindrical well around the DNA due to the solvent mediated mean-field potential or specific 
short-range ion-DNA interactions, can also replace the specific "bonding" of ions to the 
DNA surface. 

Other effects not accounted for in the dielectric continuum model are the dielectric 
discontinuity and dielectric saturation effects. The former effect emerges due to the po- 
larization of the DNA surface and affects the ion distribution outside the DNA core [183]. 
It rapidly drops off for large distances from the surface [37]. The latter effect is related 
to the water anisotropy near the DNA surface [184] and can be accounted for through a 
distance-dependent e in the electrostatic potentials [20,77,92,185-190]. A decade ago there 
was a widely accepted perception that the attraction between the DNA molecules cannot 
be explained by the electrostatic forces [81,191,192]. However the theoretical investigations 
and simulations of strongly correlated systems do not support this claim [3,20]. Finally, we 
neglect the salt-induced decrease of the DNA rigidity in salted solutions [61,101,193-197] 
and dielectric permeability of the solvent [175]. To include all these additional complications 
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into a simulation remains a very challenging task for the future. 
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TABLES 



TABLE I. Parameter sets used for the simulations of DNA-DNA interactions. 



Set qc Qs 

1 1 1 

2 2 1 

3 3 1 

4 1 2 

5 3 2 
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FIGURES 




FIG. 1. A typical snapshot of the simulation box. The DNA molecules are drawn according 
to the MAM. Black spheres on the DNA strands represent the phosphate charges. Internal grey 
spheres between the phosphates and the DNA cylindrical core are neutral. Positive (negative) salt 
ions spreaded across the simulation volume are shown as open (hatched) spheres. 
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FIG. 3. Three typical configurations for the two parallel DNA molecules when their phosphates 
charges are close to each other. The xy-plane cross-sections of DNA molecules are shown for the 
ECM. Note that only the neighboring phosphates in the inter-DNA area are shown, see the dark 
small circles labeled by letters A, B, C, D. (a) 00 + n4>s = tt, = 7r,7r/5. (b) (/>o + n(l)s = vr, 
(j) = TT — (/)s/2,7r/5 — </>s/2. (c) <j)Q + ncps = vr it </>s/2, (j) = 7r,7r/5. Here n is an integer number, 
n = 0, ±1, ±2, .... In (b) and (c) the pair of phosphates on each cylinder pertain to the same strand 
and have different z coordinates: (b) za = zb — 1.7A, Zc = zb + 1.7A; (c) za = zd, zb = zq, 

ZA — Zc = 3.4A. 
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22 25 28 31 34 37 40 

/?(A) 

FIG. 4. Reduced DNA-DNA interaction force F/ Fq versus separation distance R for the mono- 
valent counterions and monovalent salt ions (parameter set 1 of Table 1) . The unit of the force is 
Fq = ksT/P, where P is the DNA pitch length. Different salt densities are shown: Cs=0 mol/1 
(solid line), 0.024 mol/1 (dashed hue), 0.097 mol/1 (dot-dashed hue), 0.194mol/l (sohd hue with 
symbols), 0.71 mol/1 (dashed line with symbols), (a)- CM, (b)- ECM, (c)- MAM. The inset in (c) 
shows the force-salt non-monotonicity at the separation distance i?=25A. 
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FIG. 5. (a) Reduced electrostatic force F2/FQ and (b) entropic force F^/Fq components of the 
total interaction force F{R) for the MAM, compare Fig. 4c. A similar trend is observed for the 
CM and ECM. 
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FIG. 6. Reduced DNA-DNA interaction force F/Fq versus separation distance R for divalent 
counterions and monovalent salt ions (parameter set 2 of Table I) . The notation is the same as in 
Fig. 4. The inset in (c) shows the force-salt non-monotonicity at the separation distance R=2lk.. 
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FIG. 12. Reduced DNA-DNA interaction force F/Fq versus separation distance R for a divalent 
salt and monovalent counterions (parameter set 4 of Table I). The notation is the same as in Fig. 4. 
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FIG. 13. Reduced DNA-DNA interaction force F/Fq versus separation distance R for a divalent 
salt and trivalent counterions (parameter set 5 of Table I) and for the MAM. The notation is the 
same as in Fig. 4. 
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